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1.  Introduction 


Traumatic  Brain  Injury  (TBI)  is  a  major  health  issue  that  leads  to  a  range  of  neurological  and 
psychological  impairments.  In  recent  operations,  hundreds  of  thousands  of  U.S.  military 
personnel  sustained  a  mild  form  of  TBI,  often  after  exposure  to  an  explosive  blast.  Mild  TBI  is 
hard  to  diagnose  because  it  often  happens  without  external  injuries,  and  consequently  mild  TBI 
may  not  be  identified  until  the  service  member  returns  home  (7). 

The  blast  shock  wave  of  primary  TBI  is  a  transient  supersonic  pressure  wave  with  a  rapid,  less 
than  1-ms  increase  in  pressure  (compression),  density,  and  temperature,  which  corresponds  to  the 
peak  overpressure  in  figure  1,  followed  by  a  nonlinear,  microsecond  phase  of  low  pressure  or 
tension.  Such  pressure/temperature  (p  versus  t)  curves  can  be  mathematically  described  by  the 
so-called  Friedlander  waveform  (figure  1).  The  integration  of  this  curve  gives  the  positive 
impulse  per  unit  area  /,  which  is  a  useful  measure  for  positive  overpressure  duration  (2,  3). 


Figure  1.  Mathematical  representation  of  planar  Friedlander  waveform  (5). 

In  the  majority  of  TBI  cases  the  peak  pressure  is  low.  Even  exposure  to  blasts  of  10-atm  peak 
pressure  for  a  few  milliseconds  can  result  in  death  for  unprotected  people  (4).  The  mechanisms 
initiating  TBI  as  a  direct  result  of  the  shock  wave  generated  by  an  explosion  remain  obscure. 
Although  dynamic  compression,  tension,  and  shear  stress  have  all  been  proposed  to  explain 
primary  TBI,  the  precise  mechanism  of  how  this  damage  arises  is  not  fully  understood  (5). 

While  it  is  well  known  that  exposure  of  biological  cells  to  shock  waves  causes  damage  to  the  cell 
membrane,  it  is  currently  unknown  how  cellular  structural  damage  from  realistic  blast  impact 
affects  cellular  function  at  varied  time  scales  prior  to  and  subsequent  to  brain  traumas.  It  is  still 
unclear  by  which  mechanisms  damage  is  caused  and  how  it  depends  on  physical  parameters  such 
as  shock-wave  velocity,  shock-pulse  duration,  or  shock-pulse  shape. 
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Experimentally,  it  is  very  difficult  to  observe  the  dynamics  of  membrane  rupture  due  to  small 
time  scales  and  length  scales  of  the  event.  Since  a  typical  cell  membrane  is  only  several 
nanometers  thick,  while  the  pressure  front  of  a  shock  wave  travels  at  supersonic  speeds  (for 
example,  the  speed  of  sound  in  water  is  1497  m/s),  the  time  scale  during  which  a  shock  wave 
interacts  with  a  cell  membrane  is  on  the  order  of  picoseconds  (2).  As  a  consequence,  it  may 
prove  particularly  helpful  to  elucidate  the  microscopic  structural  details  and  dynamics  with 
computer  simulations. 

Cell  membranes  are  mainly  composed  of  phospholipids,  each  made  of  a  hydrophilic  segment 
(“head”)  covalently  bonded  to  a  hydrophobic  segment  (“tail”).  The  phospholipids  form  a  bilayer 
stabilized  by  hydrophobic  and  hydrophilic  interactions.  However,  simulation  of  a  realistic  cell 
membrane  is  a  challenging  task  due  to  the  structural  complexity  of  cell  membranes.  Real  brain 
cell  membranes  are  inhomogeneous  structures  that  include  multiple  lipids,  cholesterol,  and 
numerous  proteins.  The  purpose  of  this  study  is  not  to  simulate  cell  membranes  in  realistic  detail 
but  to  focus  on  a  simple  model  of  the  major  constituent  of  a  cell  membrane,  the  phospholipid 
bilayer.  In  this  work,  we  studied  the  mechanism  of  shock-wave-induced  damage  to  phospholipid 
vesicles  (liposomes)  as  the  simplified  generic  model  for  the  cell  membranes.  Vesicles  or 
liposomes  are  composed  of  phospholipids  that  encapsulate  a  core  region  of  water  inside  a 
hydrophobic  membrane.  Simulation  studies  of  vesicles  have  some  clear  advantages  over  those  of 
planar  membranes:  the  absence  of  a  periodicity  effect,  the  use  of  curvature  as  an  additional 
parameter,  and  the  closer  correspondence  to  experimental  studies  performed  with  vesicular 
systems  (6). 

The  standard  technique  to  numerically  simulate  phospholipid  bilayers  is  molecular  dynamics 
(MD).  However,  the  complexity  of  the  atomistic  force  fields  limits  them  to  membrane  patches 
and  vesicles  some  tens  of  nanometers  in  extent  for  times  of  tens  of  microseconds  or  less.  It  is 
questionable  whether  such  small  systems  can  quantitatively  reproduce  real  damage  processes  for 
several  reasons,  including  artificial  stabilization  of  the  membrane  patch  due  to  correlation 
effects,  absence  of  membrane  defects,  and  curvature,  where  membrane  rupture  likely  originates 
from  a  less-than  ideal  flat  surface.  The  short  simulation  time  of  these  models  due  to  a  small 
simulation  box  makes  it  impossible  to  study  the  time  evolution  of  the  shocked  membranes  where 
the  shock  front  can  travel  for  a  long  period  of  time  (2). 

Various  coarse-grained  (CG)  simulation  techniques,  including  inverse  Monte  Carlo  schemes, 
force-matching  approaches,  and  calibrated  techniques  based  on  thermodynamic  data,  have  been 
developed  for  phospholipid  membranes  that  allow  them  to  reach  length  and  time  scales 
inaccessible  to  atomistic  MD  by  grouping  atoms  into  particles,  thereby  reducing  the  number  of 
particles,  and  fast  dynamics  make  simulations  two  to  three  orders  of  magnitude  more  efficient 
than  the  fully  atomistic  simulation.  Advantages  of  these  approaches  include  a  relatively  high 
level  of  accuracy  and  a  closer  resemblance  to  atomistic  simulations  (<5). 
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Dissipative  Particle  Dynamics  (DPD)  (7-10)  is  a  unique  CG  technique  for  complex  fluids  that 
allows  structural  and  dynamic  relaxation  to  occur  104  times  faster  than  a  fully  atomistic 
simulation  (10).  DPD  captures  the  mass  and  momentum  conservation  that  are  responsible  for  the 
correct  hydrodynamic  behavior  and  allows  us  to  simulate  the  kinetics  of  the  liposome  system.  On 
the  contrary,  another  fast  simulation  method,  dynamic  self-consistent  field  theory,  does  not 
include  hydrodynamic  effects  that  may  influence  the  dynamics  of  structure  evolution  in  the 
system. 

DPD  simulation  of  a  CG  model  of  the  lipid  vesicle  is  used  here  to  elucidate  the  generic 
principles  of  the  cellular  membrane  damage  induced  by  the  shock  waves. 


2.  Model  and  Methods 


The  DPD  approach,  which  is  a  mesoscale  technique  for  the  complex  fluids,  allows  simulation  of 
correct  hydrodynamic  behavior.  A  DPD  system  is  composed  of  CG  soft  particles,  each 
representing  clusters  of  molecules  rather  than  individual  atoms,  moving  continuously  in  space 
and  discreetly  in  time  according  to  Newton’s  equations  of  motion.  In  a  DPD,  particles  interact 
with  each  other  via  a  pair-wise,  two-body,  short-ranged  force,  F,  that  is  written  as  the  sum  of  a 
conservative  force,  F  ,  dissipative  force,  F  ,  and  random  force,  F  : 

r=2X+IX+IX  d) 

j*i  j*i  j*i 

In  this  report,  Fc  includes  a  soft  repulsion  force,  F  ,  acting  between  two  particles  and  a 

tt 

harmonic  force,  F  ,  acting  between  adjacent  particles  in  a  polymer  chain. 

F  forces  are  derived  from  interaction  potentials  described  in  the  following.  The  remaining  two 
forces,  F  and  F  ,  provide  the  DPD  thermostat.  The  dissipative  force  modifies  the  kinetic  energy 
of  the  particles  and  is  balanced  by  the  random  force  according  to  the  fluctuation-dissipation 
theorem.  A  detailed  description  of  the  DPD  thermostat  can  found  elsewhere  (7-10). 

DPD-reduced  units  are  adopted  for  the  convenient  expression  of  parameters:  length  in  rc,  energy 
in  kBT/rc,  mass  in  particle  mass  m,  and  time  in  r  =  r  /  kBT  .  Standard  values  in  DPD  units  of 
m  =  kBT  =  1  and  overall  particle  density  pDPD  =  3  were  chosen. 

The  systems  considered  here  are  built  up  from  three  bead  species:  lipid  head  (H),  lipid  tail  (T), 
and  water  beads  (W).  The  model  lipids  have  a  headgroup  consisting  of  three  H  beads  and  two 
hydrophobic  tails,  each  of  which  is  made  up  from  four  T  beads  (figure  2)  and  corresponds  to  a 
CG  model  dimyristoyl-phosphatidylcholine  (DMPC)  of  H3(C4)2  architecture.  Each  chain-bead  C 
represents  3.5  CH2  groups  of  the  level  of  coarse-graining  Nm  =  3.5.  This  CG  model  was  used  by 
Grafmiiller  et  al.  (11). 
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Figure  2.  Snapshot  of  the  lipid  vesicle  (liposome):  (a)  the 
liposome  is  composed  from  5551  lipids,  and  (b)  a 
generic  membrane  protein  is  inserted  in  the  membrane. 
H,  T,  and  W  particles  are  colored  red,  yellow,  and  blue, 
respectively.  The  generic  rod  and  coil  chains  are 
colored  cyan  and  green,  respectively. 


The  pair  interaction  between  topologically  nonconnected  particles  is  described  by  a  weakly 
repulsive  potential: 


U DPo{r ) 


1-- 
.  TcJ 


(2) 


where  cidpd  is  the  maximum  repulsion  and  rc  =  1  is  the  cutoff  radius,  a  dpi)  differ  for  different 
bead  species  and  their  values  are  taken  from  Grafmiiller  et  al.  (11). 

Table  1.  Interaction  parameters  a,jDPD  between  pairs  of  beads,  i  and  j  could  belong  to  head  (H), 
tail  (T),  or  water  (W)  beads,  respectively. 


Bead  Type 

H 

T 

W 

H 

30 

35 

30 

T 

35 

10 

75 

W 

30 

75 

25 

In  general,  the  simulation  parameters  are  chosen  to  match  the  mesoscopic  behavior  of  the 
system. 

Topologically  connected  beads  in  lipids  interact  according  to  the  harmonic  potential,  Unarm- 

Uharm{r)  =  ^{r-ro)2  (3) 
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where  k  is  the  harmonic  constant  and  r0  is  the  equilibrium  bond  distance;  the  values  of  k  =  225 
and  ro=  0.85  are  used. 

In  addition,  the  hydrocarbon  chains  (tails)  are  stiffened  by  a  bending  potential  for  two 
consecutive  bonds  of  the  form 

UbenAr)  =  ke(l-cosd)  (4) 


where  the  potential  parameters  are  ke  =  2  and  0  =  1 80 . 

Initial  configurations  were  generated  by  randomly  placing  water  particles  into  the  simulation 
box.  The  vesicle  structure  was  preassembled  from  5551  lipids,  putting  the  head  groups  in  the 
inner  and  outer  parts  of  the  spherical  membrane.  The  final  configurations  were  obtained  by  the 
initial  equilibration  of  hydrophobic  tails  by  freezing  of  the  head  beads  and  subsequent  short 
simulation  of  fully  mobile  lipids  (figure  2a).  The  vesicle  also  contains  10  generic  membrane 
proteins  modeled  by  the  rod-coil  chains  (figure  2b).  The  beads  that  made  up  the  rigid  part  are 
hydrophobic,  and  their  stiffness  is  controlled  by  a  harmonic  angular  potential, 
k 

U a\r)~-^(fP~(Po)  >  where  %  =180  and  kip  =  50  .  The  hydrophilic  beads  belong  to  the  flexible 
part  of  the  protein. 

In  this  study,  all  of  the  equilibration  stages  were  performed  with  T  =  1.0.  The  time  step  was 

AtDPD  =  0.01. 

Shock  waves  were  initiated  with  a  momentum  reflecting  mirror,  which  is  a  standard  method  for 
generating  shock  waves  in  computer  simulation  studies  of  lipid  bilayers  (11,  12).  This  method 
allows  for  generation  of  the  fast  shock  waves  with  a  well-defined  shock-wave  front  and  good 
numerical  stability.  Infinite-mass  piston  moves  in  z  positive  direction  and  creates  shock  waves. 
All  particles  coming  into  contact  with  the  piston  surface  are  reflected,  and  the  velocity  relative  to 
the  moving  wall  is  flipped  in  z  (2). 

Upon  impact,  the  target  material  is  compressed  and  the  resulting  steep  density  gradient  initiates  a 
shock  wave.  After  a  short  initial  time,  the  piston  was  stopped  while  the  initiated  shock  wave 
continued  to  travel  further  along  z  direction  (figure  3).  The  speed  of  the  piston,  vp,  and  the  final 
piston  position  were  varied  to  get  different  shock  wave  profile.  The  pressure  was  measured  as 
function  of  time  and  distance.  The  liposome  was  placed  far  enough  in  the  simulation  box  to 
ensure  that  it  is  not  to  be  hit  by  the  piston.  This  method  is  similar  to  standard  shock  wave 
experiments  where  a  static  target  material  is  hit  by  a  fast-moving  impactor  (5). 
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Figure  3.  Snapshot  of  the  simulation  box  after  initial  compression. 

To  compare  our  simulation  results  with  experimental  data,  it  is  necessary  to  define  energy, 
length,  and  time.  The  basic  length  scale,  i.e.,  a  physical  size  of  interaction  radius,  r  =0.69nm,  is 

computed  from  r.  =  \jpDPDNmvw  nm,  where  vm  is  the  volume  of  one  water  molecule  (10,  11).  The 

coupling  of  the  time  scales  was  determined  from  a  “speed  of  sound”  in  the  DPD  fluid,  cdpd, 
which  was  found  to  be  in  a  range  of  3.818  to  4.298  DPD  units  (13).  Having  cwater  =  1497  m/s,  the 

r  c 

time  scale  was  found  to  be  equal  to  z  =  c  wa,er  =  0.25  ms.  The  same  technique  is  used  for 

CDPD 

mapping  pressure  in  this  simulation  to  the  physical  units,  where  the  equilibrium  pressure  of 
22.75  in  DPD  units  was  found  to  be  approximately  equal  to  an  atmospheric  pressure  of  14.7  psi. 

To  determine  the  speed  of  the  shock  wave,  the  z  location  of  the  shock  wave  front  was  identified 
as  a  function  of  simulation  time.  The  position  of  the  liposome  was  monitored  during  the 
simulation,  and  pressure  impulse  and  profile  were  computed. 

Since  study  of  the  time  evolution  of  the  shocked  vesicle  requires  a  large  box  along  the  direction 
of  the  shock  impulse,  Lz,  to  allow  the  shock  front  to  travel  for  a  long  period  of  time,  the 
simulation  box  is  relatively  large — 35  x  35  x25  nm — and  the  liposome  is  20  nm  in  diameter.  The 
simulation  box  contains  approximately  2,500,000  particles. 

All  simulations  were  executed  using  LAMMPS  software  from  Sandia  National  Laboratories  (14, 
15). 


3.  Results  and  Discussion 


3.1  Simulation  of  the  Blast  Waves  with  Low  Peak  Pressure 

Figure  4  shows  an  example  of  pressure  profiles  during  the  simulation  at  various  distances  from 
the  initiation  of  the  shock  wave  that  mimics  the  distance  of  the  center  of  the  explosion.  The  point 
of  the  shock  wave  initiation  is  defined  here  as  the  center  of  the  explosion.  As  expected,  the  peak 
overpressure  decreases  with  distance  from  the  center  of  the  explosion  (see  figure  4a).  These 
curves  closely  mimic  the  supersonic  shock  wave  produced  by  high  explosives  passing  through 
the  head.  This  wave  causes  almost  instantaneous  overpressure,  followed  by  a  longer  wave  of 
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underpressure,  and  a  second  wave  of  overpressure  (/).  The  similar  shock  wave  was 
experimentally  obtained  by  a  pneumatic  device  that  delivers  shock  waves  (5). 

The  curves  shown  in  figure  4  correspond  to  the  shock  wave  with  the  Mach  number  M  =  1.29, 
where  the  Mach  number  is  a  dimensionless  quantity  representing  the  ratio  of  speed  of  the  shock 
wave  in  a  fluid  and  the  local  speed  of  sound.  The  magenta  curve  represents  the  change  in 
pressure  at  the  location  of  the  liposome,  which  was  placed  at  1 17.4  nm  from  the  center  of  the 
explosion  and  did  not  move  significantly  (figure  4a).  The  pressure  profiles  that  corresponds  to 
the  same  system  at  the  time  of  the  peak  overpressure  at  the  vesicle  location  (t  =  6.75  ms)  is 
shown  in  figure  4b. 


Figure  4.  (a)  The  supersonic  shock  blast  wave  at  the  various 
distances  from  its  launch.  The  liposome  location 
corresponds  to  the  magenta  curve.  The  second  peak  in 
blue  and  magenta  curves  represents  the  second  positive 
phase  due  to  reflection.  The  Mach  number  is  1.29. 

(b)  The  pressure  profile  at  time  =  6.75  ms  after  the 
explosion  at  the  liposome  location. 
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To  get  different  pressure  profiles  for  mild  shock,  the  piston  speed  vp  was  varied,  keeping  the 
compression  of  the  box  constant.  Values  of  vp  were  equal  to  7485  m/s,  14,970  m/s,  and 
22,455  m/s,  respectively.  This  technique  produces  the  blast  waves  with  the  Mach  numbers  of 
1.29,  1.43,  1.48,  respectively,  and  overpressure  peaks  equal  to  17.8,  24.6,  and  31.4  psi, 
respectively.  The  overpressure  peaks  were  measured  at  approximately  117.6  nm  (at  the  location 
of  the  liposome). 

The  results  are  reported  in  figure  5.  The  liposome  is  moved  in  the  positive  z  direction  during  the 
simulation,  as  the  shock  wave  hits  it,  and  then  the  liposome  flows  in  the  negative  z  direction  as 
the  reflective  wave  reaches  its  location  (the  second  positive  phase).  The  movement  of  liposome 
in  the  lateral  direction  to  the  shock  wave  is  found  to  be  negligible  (less  than  1%).  With  the 
increase  of  vp,  the  reflective  wave  hits  the  liposome  before  the  pressure  reaches  the  ambient 
value. 


time,  ms 


Figure  5.  (a)  Shock  wave  at  the  location  of  the  liposome  for  the 
various  piston  speeds  vp  (b)  Movement  of  the  center  of 
mass  of  the  liposome  in  z  direction. 
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To  characterize  the  deviation  from  an  equilibrium  vesicle  structure,  the  diagonal  principle 
moments  of  the  gyration  tensor,  Xxx ,  ,  and  X:, ,  were  computed,  which  are  given  by 

4  =  T7  Sfe  "  Vcom  )(rid  ~  rcom  )  (5) 

iV  i= 0 

where  rcom  represents  the  coordinates  of  the  liposome  center  of  mass  and  N  is  number  of  beads  in 
the  liposome.  Figure  6  shows  the  change  in  in  time.  The  results  of  this  study  suggest  that 

changes  in  Xxx  and  X2yy  are  insignificant  compared  with  the  change  in  .  The  liposome  is 
compressed  during  the  positive  phase,  then  it  is  extended  at  the  negative  phase  and  the  process 
repeats  with  lower  values  of  compression/extension  for  the  reflective  wave.  Notably,  an  extent  of 
the  liposome  compression  does  not  strongly  depend  on  the  pressure  pulse  compared  with  an 
extent  of  liposome  extension  (figure  6b).  This  finding  suggests  that  the  liposome  extension  is 
very  sensitive  to  the  change  in  the  negative  pressure  phase,  which  can  consequently  lead  to  more 
structural  damage  of  the  cell  membranes. 
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Figure  6.  Change  in  the  diagonal  principle  moments  of  the  gyration 
tensor  (a)  )}xx ,  ,  and  (b)  X~_  during  the  simulation. 

Figure  7a  shows  vesicle  structures  at  the  peak  overpressure,  which  corresponds  to  the  most 
compressed  liposome  configurations.  At  this  range  of  overpressure  (17.8-31.4  psi),  the  liposome 
membrane  does  not  sustain  any  significant  damage;  even  the  liposome  is  compressed.  Figure  7b 
demonstrates  substantial  extension  of  the  vesicle  at  the  considering  negative  pressure  range 
(1. 7-6.0  psi). 
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Figure  7.  Snapshot  from  the  liposome  simulation  at  the  various 

piston  speeds:  (a)  compressed  structures  and  (b)  extended 
structures.  For  visual  clarity,  water  particles  are  not  shown. 

H  and  T  particles  are  colored  red  and  yellow,  respectively. 

The  generic  rod  and  coil  chains  are  colored  cyan  and  green, 
respectively. 

Even  visually  intact,  the  liposome  membrane  starts  to  have  perforation  in  its  hydrophobic  (inner) 
layer  at  negative  pressure,  approximately  -4.8  psi.  Figure  8  shows  a  small  hole  in  hydrophobic 
part  of  the  bilayer.  The  results  imply  that  the  cellular  membrane  permeability  is  altered  even  for 
the  small  negative  pressure,  and  consequently  sodium,  potassium,  and  calcium  ion  influx  may 
change.  This  finding  suggests  that  the  cellular  damage  leading  to  the  TBI  may  take  place  because 
of  the  relatively  low  negative  pressure. 


Figure  8.  Part  of  the  liposome 

hydrophobic  layer  with  a 
hole.  For  visual  clarity,  only 
hydrophobic  particles  are 
shown. 

3.2  Simulation  of  the  Blast  Waves  with  High  Peak  Pressure 

In  order  to  get  higher  overpressure,  two  simulations  were  performed  with  the  piston  speed  vp 
equal  to  7485  m/s  and  14,970  m/s,  respectively.  In  these  simulations,  the  piston  was  allowed  to 
move  longer,  and  therefore  both  positive  and  negative  pressure  phase  values  are  greater  than  the 
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previous  simulation  of  the  mild  shock  (figures  5  a  and  9a).  The  Mach  numbers  of  these  waves  are 
close:  1.85  and  1.89,  which  correspond  to  vp  equal  to  7485  m/s  and  14,970  m/s,  respectively. 
Since  the  shock  waves  move  faster  in  this  series  of  simulations,  the  reflection  wave  reaches  the 
liposome  earlier  than  in  the  previous  simulations.  The  durations  of  the  positive  and  negative 
phases  are  close  to  the  previous  simulation  (figures  5a  and  9a). 


The  liposome  flows  in  the  positive  z  direction  at  the  greater  distance  compared  with  the  mild 
shock  (figures  5b  and  9b). 
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Figure  9.  (a)  Shock  wave  at  the  location  of  the  liposome  for  the 
various  piston  speeds  vp.  (b)  Movement  of  the  center  of 
mass  of  the  liposome  in  z  direction. 

Change  in  the  diagonal  principle  moments  of  the  gyration,  Xxx  and  ,  is  slightly  more 
pronounced  than  the  mild  shock.  Xzz  decreases  significantly  more  for  this  series  due  to  higher 
peak  overpressure  (figures  6b  and  10b).  Interestingly,  that  the  liposome  undergoes  less  extension 
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in  the  z  direction  for  the  simulation  with  vp  =  14,970  m/s  than  the  liposome  in  the  simulation  of 
vp  =  7485  m/s  (figure  10b).  This  effect  caused  by  the  reflective  wave  that  reaches  the  liposome 
before  the  vesicle  was  fully  stretched. 


Figure  10.  Change  in  the  diagonal  principle  moments  of  the 

gyration  tensor  (a)  ,and(b)  7,7 during  the 

simulation. 

The  simulation  results  show  that  perforation  in  the  hydrophobic  (inner)  layer  of  the  liposome 
membrane  begins  at  the  positive  pressure.  Note  that  approximately  90  psi  of  the  overpressure 
causes  the  significant  membrane  perforation.  Figure  1 1  shows  a  large  hole  in  hydrophobic  layer 
of  the  compressed  vesicle. 
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Figure  1 1 .  Compressed  liposome.  H  and  T  particles  are  colored  red 
and  yellow,  respectively.  For  visual  clarity,  water  and 
protein  particles  are  not  shown.  The  circle  shows  a 
location  of  a  hole  in  the  hydrophobic  layer. 


4.  Conclusion 


The  simulation  findings  suggest  that  the  cellular  membrane  permeability  is  altered  even  for  a 
relatively  mild  shock  wave  (M  <  1.5)  in  the  negative  pressure  phase.  It  was  found  that  the  small 
negative  pressure,  approximately  -5  psi,  may  lead  to  the  cellular  damage  by  changing  the  ion 
influx.  The  faster-moving  shock  waves  (M  ~  1.9),  which  gave  higher  overpressure  peaks,  already 
cause  the  cellular  damage  in  the  positive  pressure  phase  if  peak  overpressure  is  greater  than  60 
psi.  In  future  work,  we  will  extend  our  approach  to  characterize  the  change  of  specific 
permeability  of  the  cellular  membranes  induced  by  the  shock  blast  waves.  Since  the  real  brain 
cell  membranes  are  inhomogeneous  structures  that  include  multiple  lipids,  cholesterol,  and 
numerous  proteins,  we  will  develop  a  model  of  the  cellular  membrane  with  a  higher  level  of 
accuracy.  The  results  of  this  improved  model  of  the  cellular  membrane  will  be  incorporated  into 
the  multilevel  modeling  of  electrochemical  connectivity  of  large-scale  brain  networks. 
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